Peak heatmap
Bing Shui
2021-02-02
Prepare libraies
suppressMessages(
suppressWarnings(
c(library(data.table),
library(tidyverse),
library(rtracklayer),
library(DESeq2),
library(Biostrings),
library(genomation),
library(RColorBrewer),
library(CLIPanalyze))
)
)## [1] "data.table" "stats" "graphics"
## [4] "grDevices" "utils" "datasets"
## [7] "methods" "base" "forcats"
## [10] "stringr" "dplyr" "purrr"
## [13] "readr" "tidyr" "tibble"
## [16] "ggplot2" "tidyverse" "data.table"
## [19] "stats" "graphics" "grDevices"
## [22] "utils" "datasets" "methods"
## [25] "base" "rtracklayer" "GenomicRanges"
## [28] "GenomeInfoDb" "IRanges" "S4Vectors"
## [31] "BiocGenerics" "parallel" "stats4"
## [34] "forcats" "stringr" "dplyr"
## [37] "purrr" "readr" "tidyr"
## [40] "tibble" "ggplot2" "tidyverse"
## [43] "data.table" "stats" "graphics"
## [46] "grDevices" "utils" "datasets"
## [49] "methods" "base" "DESeq2"
## [52] "SummarizedExperiment" "DelayedArray" "matrixStats"
## [55] "Biobase" "rtracklayer" "GenomicRanges"
## [58] "GenomeInfoDb" "IRanges" "S4Vectors"
## [61] "BiocGenerics" "parallel" "stats4"
## [64] "forcats" "stringr" "dplyr"
## [67] "purrr" "readr" "tidyr"
## [70] "tibble" "ggplot2" "tidyverse"
## [73] "data.table" "stats" "graphics"
## [76] "grDevices" "utils" "datasets"
## [79] "methods" "base" "Biostrings"
## [82] "XVector" "DESeq2" "SummarizedExperiment"
## [85] "DelayedArray" "matrixStats" "Biobase"
## [88] "rtracklayer" "GenomicRanges" "GenomeInfoDb"
## [91] "IRanges" "S4Vectors" "BiocGenerics"
## [94] "parallel" "stats4" "forcats"
## [97] "stringr" "dplyr" "purrr"
## [100] "readr" "tidyr" "tibble"
## [103] "ggplot2" "tidyverse" "data.table"
## [106] "stats" "graphics" "grDevices"
## [109] "utils" "datasets" "methods"
## [112] "base" "genomation" "grid"
## [115] "Biostrings" "XVector" "DESeq2"
## [118] "SummarizedExperiment" "DelayedArray" "matrixStats"
## [121] "Biobase" "rtracklayer" "GenomicRanges"
## [124] "GenomeInfoDb" "IRanges" "S4Vectors"
## [127] "BiocGenerics" "parallel" "stats4"
## [130] "forcats" "stringr" "dplyr"
## [133] "purrr" "readr" "tidyr"
## [136] "tibble" "ggplot2" "tidyverse"
## [139] "data.table" "stats" "graphics"
## [142] "grDevices" "utils" "datasets"
## [145] "methods" "base" "RColorBrewer"
## [148] "genomation" "grid" "Biostrings"
## [151] "XVector" "DESeq2" "SummarizedExperiment"
## [154] "DelayedArray" "matrixStats" "Biobase"
## [157] "rtracklayer" "GenomicRanges" "GenomeInfoDb"
## [160] "IRanges" "S4Vectors" "BiocGenerics"
## [163] "parallel" "stats4" "forcats"
## [166] "stringr" "dplyr" "purrr"
## [169] "readr" "tidyr" "tibble"
## [172] "ggplot2" "tidyverse" "data.table"
## [175] "stats" "graphics" "grDevices"
## [178] "utils" "datasets" "methods"
## [181] "base" "CLIPanalyze" "GenomicAlignments"
## [184] "Rsamtools" "GenomicFeatures" "AnnotationDbi"
## [187] "plyr" "RColorBrewer" "genomation"
## [190] "grid" "Biostrings" "XVector"
## [193] "DESeq2" "SummarizedExperiment" "DelayedArray"
## [196] "matrixStats" "Biobase" "rtracklayer"
## [199] "GenomicRanges" "GenomeInfoDb" "IRanges"
## [202] "S4Vectors" "BiocGenerics" "parallel"
## [205] "stats4" "forcats" "stringr"
## [208] "dplyr" "purrr" "readr"
## [211] "tidyr" "tibble" "ggplot2"
## [214] "tidyverse" "data.table" "stats"
## [217] "graphics" "grDevices" "utils"
## [220] "datasets" "methods" "base"Load miR-peaks dataset
mirs.peaks <- readRDS("Datafiles/miRNA-peaks-list-12232019.rds")Analysis
Prepare bw file for visualization
# combine bigwig files of the same genotype
$ bigWigMerge Exp-HF1_LIB044916_GEN00173026_R1_barcode.uniq.bw Exp-HF2_LIB044916_GEN00173029_R1_barcode.uniq.bw Exp-HF3_LIB044916_GEN00173032_R1_barcode.uniq.bw HF_merged.bedGraph
$ bigWigMerge Exp-HFK1_LIB044916_GEN00173035_R1_barcode.uniq.bw Exp-HFK2_LIB044916_GEN00173038_R1_barcode.uniq.bw Exp-HFK3_LIB044916_GEN00173041_R1_barcode.uniq.bw HFK_merged.bedGraphconvert chromosome name from Ensembl to UCSC format.
chromosomename <- read.table("../Joint_analysis/Ensembl.mm10.chrom.sizes.txt", sep = "\t")[,1][1:22]
hf_bed <- read.table("../Joint_analysis/HF_merged.bedGraph", sep = "\t")
hf_bed <- hf_bed %>% filter(V1 %in% chromosomename)
hf_bed$V1 <- paste0("chr", hf_bed$V1)
hf_bed$V1[hf_bed$V1 == "chrMT"] <- "chrM"
write_delim(hf_bed, file = "../Joint_analysis/HF_merged_edited.bedGraph", delim = "\t", col_names = F)
hfk_bed <- read.table("../Joint_analysis/HFK_merged.bedGraph", sep = "\t")
hfk_bed <- hfk_bed %>% filter(V1 %in% chromosomename)
hfk_bed$V1 <- paste0("chr", hfk_bed$V1)
hfk_bed$V1[hfk_bed$V1 == "chrMT"] <- "chrM"
write_delim(hfk_bed, file = "../Joint_analysis/HFK_merged_edited.bedGraph", delim = "\t", col_names = F)# bedGraph to bigWig conversion
$ awk 'NR!=1' HF_merged_edited.bedGraph > input.HF_merged.bedGraph
$ sort -k1,1 -k2,2n input.HF_merged.bedGraph > sorted.input.HF_merged.bedGraph
$ bedGraphToBigWig sorted.input.HF_merged.bedGraph mm10.chrom.sizes HF_merged.bw
$ awk 'NR!=1' HFK_merged_edited.bedGraph > input.HFK_merged.bedGraph
$ sort -k1,1 -k2,2n input.HFK_merged.bedGraph > sorted.input.HFK_merged.bedGraph
$ bedGraphToBigWig sorted.input.HFK_merged.bedGraph mm10.chrom.sizes HFK_merged.bwpeaks <- mirs.peaks
peaks <- lapply(peaks, FUN = function(x) {x[order(x$count, decreasing = T)]})
peaks <- lapply(peaks, FUN = function(x) {x[x$padj<=0.01]})
mybw.dir <- "../Joint_analysis"
mybw.files <- list.files(mybw.dir, pattern = "bw$", full.names = T)[7:8]
print(mybw.files)## [1] "../Joint_analysis/HF_merged.bw" "../Joint_analysis/HFK_merged.bw"cols <- brewer.pal(name = "Set2", n = 8)
reds <- brewer.pal(name = "Reds", n = 9)
blues <- brewer.pal(name = "Blues", n = 9)
mycolors <-c(blues[6], reds[6]) #HF, HFK
light.colors <- alpha(c(blues[3], reds[2]), 0.4)Histogram
histogram plot function
#plot histograms
peaks_meta <- function(mypeaks = peaks,
miRNA_family = "miR-451a",
dispersion = "se",
dispersion.col = NULL,
coordinates = c(-400, 400),
line.col = mycolors,
winsorize = c(0,99),
title = ""){
suppressWarnings(mypeaks <- lapply(mypeaks, FUN = function(x) {resize(x, width = sum(abs(coordinates)), fix="center")}))
mypeaks <- GRangesList(mypeaks)
mysml <- ScoreMatrixList(targets=mybw.files, window=mypeaks[[miRNA_family]], type = "bigWig", strand.aware = T)
mysampleInfo <- data.frame(basename(mybw.files), c("HF", "HFK"))
names(mysampleInfo) = c("sample", "genotype")
names(mysml) = mysampleInfo$genotype[match(names(mysml), mysampleInfo$sample)]
plotMeta(mysml, profile.names = c("HF", "HFK"), xcoords = coordinates, dispersion = dispersion, main =title, line.col = line.col, winsorize = winsorize, dispersion.col = dispersion.col)
}len <- sapply(mirs.peaks, function(x) length(x))
mirs.peaks <- mirs.peaks[order(-len)]
mirna <- names(mirs.peaks)
for (i in 1:10) {
peaks_meta(mypeaks = peaks, dispersion = NULL, miRNA_family = mirna[i], title = mirna[i])
peaks_meta(mypeaks = peaks, dispersion = "se", miRNA_family = mirna[i], title = mirna[i],
dispersion.col = light.colors)
}
Heatmap
plot heatmaps
peaks_heat <- function(mypeaks = peaks,
miRNA_family = "miR-451",
col = blues9,
coordinates = c(-400, 400),
order_rows = F,
winsorize_parameters = c(1,98)){
suppressWarnings(mypeaks <- lapply(mypeaks, FUN = function(x) {resize(x, width = sum(abs(coordinates)), fix="center")}))
mypeaks <- GRangesList(mypeaks)
mysml <- ScoreMatrixList(targets=mybw.files[c(1,2)], window=mypeaks[[miRNA_family]], type = "bigWig", strand.aware = T)
mysampleInfo <- data.frame(basename(mybw.files), c("HF", "HFK"))
names(mysampleInfo) = c("sample", "genotype")
names(mysml) = paste(mysampleInfo$genotype[match(names(mysml), mysampleInfo$sample)], miRNA_family , sep = "_")
mysml.scaled = scaleScoreMatrixList(mysml)
#multiHeatMatrix(mysml.scaled, xcoords = coordinates, col = col)
multiHeatMatrix(mysml, common.scale = T, xcoords = coordinates, winsorize = winsorize_parameters, col = col, order = order_rows)
}colfunc <- colorRampPalette(c("white", "blue"))
mycols <- colfunc(128)for (i in 1:10) {
peaks_heat(mypeaks = peaks, miRNA_family = mirna[i], col = blues9, winsorize_parameters = c(0,99))
}
SessionInfo
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.4
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] CLIPanalyze_0.0.10 GenomicAlignments_1.24.0
## [3] Rsamtools_2.4.0 GenomicFeatures_1.40.1
## [5] AnnotationDbi_1.50.3 plyr_1.8.6
## [7] RColorBrewer_1.1-2 genomation_1.20.0
## [9] Biostrings_2.56.0 XVector_0.28.0
## [11] DESeq2_1.28.1 SummarizedExperiment_1.18.2
## [13] DelayedArray_0.14.1 matrixStats_0.58.0
## [15] Biobase_2.48.0 rtracklayer_1.48.0
## [17] GenomicRanges_1.40.0 GenomeInfoDb_1.24.2
## [19] IRanges_2.22.2 S4Vectors_0.26.1
## [21] BiocGenerics_0.34.0 forcats_0.5.1
## [23] stringr_1.4.0 dplyr_1.0.4
## [25] purrr_0.3.4 readr_1.4.0
## [27] tidyr_1.1.2 tibble_3.0.6
## [29] ggplot2_3.3.3 tidyverse_1.3.0
## [31] data.table_1.13.6
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-0 ellipsis_0.3.1 fs_1.5.0
## [4] rstudioapi_0.13 farver_2.0.3 bit64_4.0.5
## [7] lubridate_1.7.9.2 xml2_1.3.2 splines_4.0.3
## [10] cachem_1.0.1 impute_1.62.0 geneplotter_1.66.0
## [13] knitr_1.31 jsonlite_1.7.2 seqPattern_1.20.0
## [16] broom_0.7.4 gridBase_0.4-7 annotate_1.66.0
## [19] dbplyr_2.0.0 compiler_4.0.3 httr_1.4.2
## [22] biosignals_0.1.0 backports_1.2.1 assertthat_0.2.1
## [25] Matrix_1.3-2 fastmap_1.1.0 cli_2.3.0
## [28] htmltools_0.5.1.1 prettyunits_1.1.1 tools_4.0.3
## [31] gtable_0.3.0 glue_1.4.2 GenomeInfoDbData_1.2.3
## [34] reshape2_1.4.4 rappdirs_0.3.3 Rcpp_1.0.6
## [37] cellranger_1.1.0 vctrs_0.3.6 xfun_0.20
## [40] rvest_0.3.6 lifecycle_0.2.0 XML_3.99-0.5
## [43] zlibbioc_1.34.0 scales_1.1.1 BSgenome_1.56.0
## [46] hms_1.0.0 curl_4.3 yaml_2.2.1
## [49] memoise_2.0.0 biomaRt_2.44.4 stringi_1.5.3
## [52] RSQLite_2.2.3 highr_0.8 genefilter_1.70.0
## [55] plotrix_3.8-1 BiocParallel_1.22.0 rlang_0.4.10
## [58] pkgconfig_2.0.3 bitops_1.0-6 evaluate_0.14
## [61] lattice_0.20-41 bit_4.0.4 tidyselect_1.1.0
## [64] magrittr_2.0.1 R6_2.5.0 generics_0.1.0
## [67] DBI_1.1.1 pillar_1.4.7 haven_2.3.1
## [70] withr_2.4.1 survival_3.2-7 RCurl_1.98-1.2
## [73] modelr_0.1.8 crayon_1.4.0 KernSmooth_2.23-18
## [76] BiocFileCache_1.12.1 rmarkdown_2.6 progress_1.2.2
## [79] locfit_1.5-9.4 readxl_1.3.1 blob_1.2.1
## [82] reprex_1.0.0 digest_0.6.27 xtable_1.8-4
## [85] openssl_1.4.3 munsell_0.5.0 askpass_1.1